A moment based approach to the dynamical solution of the Kuramoto model 
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We examine the dynamics of the Kuramoto model with a new analytical approach. By defining 
an appropriate set of moments the dynamical equations can be exactly closed. We discuss some 
applications of the formalism like the existence of an effective Hamiltonian for the dynamics. We 
also show how this approach can be used to numerically investigate the dynamical behavior of the 
model without finite size effects. 



The study of the dynamical behavior of systems with a 
very large number of mutual interacting units is a longly 
debated subject. It is a topic with interest in many dif- 
ferent interdisciplinar fields. The cooperation between 
the members of a population may lead to very rich dy- 
namical situations ranging from chaos, periodicity, phase 
locking, synchronization to self-organized critical states, 
just to cite a few In the presence of disorder such 

interaction can be frustrated and this yields new types 
of behavior. In the realm of disordered systems much 
work has been devoted to the study of models with relax- 
ational dynamics, for instance spin-glass models ||. In 
those cases there exists an Hamiltonian function which 
governs the dynamics of the system. A large body of in- 
formation can be obtained by using the tools of statistical 
mechanics. One of the main results at equilibrium is that 
the fluctuation-dissipation theorem is obeyed. But it is 
definitely interesting to study the dynamical behavior of 
dissipative systems in the presence of external driving 
forces. 

A simple model of this type was proposed by Kuramoto 
to analyze synchronization phenomena in populations of 
weakly nonlinearly coupled oscillators Q . It has become 
the subject of extensive studies in recent years due to its 
applications to biology, chemistry and physics ||. The 
purpose of this letter is to present a new analytical ap- 
proach to the Kuramoto model based on the definition 
of a suitable hierarchy of moments. It allows to repro- 
duce previous known results and, in addition, gives a 
new insight into the nature of the problem. Here, we will 
present the method and consider its potential applica- 
tions leaving detailed analysis for future work. Through 
this formalism it is possible to analyze some aspects of 
the model that deserve special attention. As an example, 
it has been suggested that, under certain conditions, it is 
possible to define a suitable Hamiltonian function from 
which it is possible to compute stationary properties of 
the system within the usual thermodynamic formalism 



such as ground states and universality classes at zero 
temperature || and equilibrium Boltzmann distribution 
in the more general case at finite temperature {jj. Our 
method can answer this question in a simple way. We will 
also show how our approach can be used to numerically 
investigate the behavior of the Kuramoto model free of 
finite size effects. Particular results will be obtained for 
the bimodal distribution case. 

Our formalism complements other recent theories de- 
veloped to analyze the Kuramoto model. In particular, 
it is worthwhile to mention the order function approach 
Q useful for studying properties of the stationary states 
of the system as well as the critical exponent of the order 
parameter at the onset of entrainment. Another interest- 
ing method was proposed in || based on kinetic theory 
and suitable to deal with questions related to the time 
dependence of the probability density of the system. 

The Kuramoto model is defined by a set of N oscilla- 
tors whose state can be specified in terms of only one de- 
gree of freedom, the phase. Each phase {(pi; 1 < i < N} 
follows the dynamical equation 
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where u)i is the intrinsic frequency of the oscillator ran- 
domly chosen from a distribution of density g(w), Kq 
is the strength of the coupling which, as in the original 
case, we will consider ferromagnetic although more com- 
plex situations have been analyzed in the literature [ p)| . 
Finally, rji (t) denotes a gaussian independent white-noise 
process 



<Vi(t)Vj{t') >=2TS ij 5(t-t'). 



(2) 



Without any other ingredient there is a competition be- 
tween the coupling, which tends to synchronize all the 
oscillators, and the noise (frequencies plus thermal noise) 
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which breaks the coherence. For a critical K c there is a 
spontaneous transition from incoherence to a new state 
where a macroscopic number of units are synchronized. 
To solve the dynamics of the Kuramoto model we define 
the the following set of moments, 



1 



AT 



< exp(ik(j>j) > uij 1 



(3) 



where i stands for the imaginary unit and k,m are in- 
tegers in the range (—00, 00), [0, 00) respectively. The 
averages < (•) > and (•) indicate averages over the noise 
and frequencies distribution respectively. The definition 
of this set of moments is the basis of the new dynamical 
approach we are proposing. Note that H™ is the more 
natural object we can construct which is invariant un- 
der the local transformation fa — > fa + 2tt. This is also 
the local symmetry of the dynamical equations (0). It is 
possible to show that the moments H™ are self-averaging 
with respect to the thermal noise 0. The equation of mo- 
tion for H™ can be easily derived using the eq.(||), 



dH" k 
dt 



K Q k 



k 2 TR m + ikH m+l 

(4) 



where we have defined the h k — H®. The term /c 2 Ti7™ 
can be simply obtained using the Gaussian representation 
for the noise rji, doing an integration by parts and using 
the regularisation condition §^py = 1/2. In this form the 

equations are closed because the time operator acting 
on the moment H ™ generates always new moments of the 
same type. If wc define the time dependent generating 
function g t (x,y) = E^-ooXm=o ex PHM|^ 
it is easy to check from eq.(0) that it satisfies the fol- 
lowing differential equation, 



definition of the generating function it is straightforward 
to check that gt(x, y) = jf S^Li ^{"Pj ~ x ) ex P(2/ w j)- F° r 
y — this is nothing else than the probability density of 
one oscillator to have phase x. In this way we recover the 
results obtained by Bonilla [[TTJ using the path integral 
formalism. The hierarchy of equations (|4|) only depends 
on the time evolution of the moment h\(h-\ = h\) which 
is the order parameter of the problem. The full set of 
moments are self-consistently computed using the con- 
ditions hi(t) — J 27r exp(ix)g t (x, 0)dx and H™ = uj m . 
In this way, we have reached a dynamical solution of the 
problem identical to that found in some mean-field glassy 
models where dynamical equations can be exactly closed 

Once we have presented the formulation we present its 
applications. The method furnishes a transparent way to 
show why a static description based on conventional equi- 
librium statistical mechanics cannot give reliable infor- 
mation about the long-time properties of the Kuramoto 
model (like the existence of stationary states). Let us 
stress from the beginning that our approach deals with 
the dynamics of the model once the thermodynamic limit 
N — > 00 is taken first than the limit t — > 00. Our results 
are meaningful in this case. In the other dynamical ap- 
proach one does the thermodynamic limit N — > 00 after 
solving the dynamics in the long-time limit t — > 00. This 
yields quite different results. These considerations in- 
terest specially to the results reported in || where the 
last situation has been analyzed. Which of the two ap- 
proaches is valid relies on the timescales one is able to ob- 
serve. Because the crossover time which matches the two 
approaches grows extremely fast with the system size, 
our approach (which , on the other hand, is the most con- 
ventional one) is more suited to investigate what can be 
observed in macroscopic systems in realistic timescales. 
Let us consider the following Hamiltonian function, 



-m=-dx-( A{x > t)9 < 



T 



d 2 gt d 2 g t 



dx 2 dxdy 



(5) 



where 
A(x,t) 



— ^(exp(— ix)h\ — exp(ix)h-i) = 
2i 



N 



N 



sin(x — <j>j) — K r sin(6* — x) (6) 



i=i 



and we have expressed hi = h*_ l = rexp(i9) where r is 
the parameter which measures the coherence (synchroni- 
sation) between oscillators. By substituting eq.(0) in the 



*The derivation of this result comes from the fact that a 
probability density of oscillators can be defined for the Ku- 
ramoto model, see Strogatz and Mirollo in p. 



Keff = 



Ko 

' N 



i) - Wi<t>i 



(7) 



where the phases <pi are restricted to the interval [— n, tt) 
in order Ti. e j f to be bounded from below. This Hamilto- 
nian function is the only candidate which generates the 

equations (Q) in the Langevin dynamics fa — g^ /J + 7 ?i- 

We will show that the equilibrium solutions of ^) are 
not stationary solutions of the dynamics eq.(|l|) even at 
zero temperature. To prove this result we compute the 
partition function of eq.(0) at temperature T = 4 and 
evaluate the moments HJT'(eq) in equilibrium. The com- 
putations are quite simple since the disorder in H e ff is 
only site dependent. For sake of simplicity we will con- 
sider here the case in which the disorder distribution 
is symmetric (g(to) — g(—u>)). It is easy to obtain the 
equilibrium values of the different moments. We obtain 
HJ! l (eq) = FJJ l (eq) exp(ik8) where 9 is an arbitrary phase 



2 



and the FJJ l (eq) are given by 



F^{eq)^uo m R%{f3K r) 



(8) 



where the A(uS) — J dujg{ijj)A(uo) is the average over the 
frequency distribution and Rf(x) = J Jl ^ . The J£(x) 
are generalized Bessel-like functions defined by, 



Vodd as a function of the temperature in the bimodal dis- 
tribution for different values of Ko This proves that 
equilibrium states of H e f / at finite temperature and also 
the ground states of TC e f / are not stationary states of the 
dynamics. Note that we cannot discard the fact that lo- 
cal minima (but not global) of H e f / at zero temperature 
are fixed points of the dynamics. 



J k( x ) = d 4> exp(i/c0 + x cos(» + j3u4>) (9) 
Jo 



where r is fixed by the condition r 



F?(eq) 



To check if the F™ (eq) are stationary solutions of the 
dynamics we plug them in eq.(^) and use the recursion 
relation, 



Jk+li x ) — Jk~li x ) " — ~Jk( x ) 



2i exp(a;) 



(exp(27r/3i 



LU ) - 1 



obtaining after some manipulations, 

dF}p(eq) 



dt 



ikv r , 



where 



TeMPK r)- 



(exp(27r f3uj) - 1 



(10) 



(ii) 



(12) 



Due to the symmetry property of the g(u)) it is easy 
to check that the v m always vanishes for m even. But it 
never does for odd m. In the high temperature regime 



/3 — ^ it is possible to show that V^m-i 



u> 



2m 



+ o(p 3 



In the limit T — > it is easy to check that V2m-i = 
uj 2m + O(^). The general result is that odd m-moments 
violate stationarity [J We conclude that the time deriva- 
tive of the equilibrium moments Hu l (eq) with odd m get 
an imaginary contribution proportional to i = exp(i : |) 
which is transverse in the complex plane to the H™(eq) 
itself. Note that the term v m is the angular velocity or 
time derivative of the global phase for all the moments 
which only depends on the number m. In the case of the 
bimodal distribution g(u) — \{5(uj — ojq) + S(uj + u>q)) 
all moments reduce to two different moments (see below) 
depending if m is even or odd. In figure 1 we show the 
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FIG. 1. Vodd as a function of T for the bimodal distribution 
for three values of Ko — 0.5, 1,2 from bottom to top. The 
cusp in v dd moves to higher temperatures as Ko increases. 

Now we want to show how equations (||) can be used 
as a powerful tool to investigate the dynamical behavior 
of the Kuramoto model. We will consider the case of 
the bimodal distribution because the phase diagram of 
the model is very rich. We have checked the situation 
for other distribution of frequencies and the results have 
been always very satisfactory. In this case the full set 
of moments H™ reduces to two different sets, fk = H™ 



for m even and gu = H™ for m odd. 
dynamical equations read in this case, 



The full set of 



dfk 
dt 



~TTT — ^—(fk+lf-1 - fk-lfl) 



k 2 Tf k 



ikg k (13) 



Ti. e f f only changes by a a constant $ cjj if all the phases 
4>i are changed to 4>i + $ 

"''An exception to this rule are the moments H™, among them 
the average frequency, because k = and therefore the right 
hand side of eq.(0) vanishes. 



dgk 
dt 



K k 



{gk+if-i - gk-xh) 



k 2 Tg k 



ikf k (14) 



§ In the rest of the paper and without loss of generality we 
will take wo = 1 
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By denning p+(x) = \ J2kL-oo> exp(-ikx)(f k + g k ) 
and p-(x) = ^J2kL-oo ex P(~ikx)(f k - g k ) we observe 
that these are the probability densities of having one os- 
cillator with phase x and natural frequencies +1 and —1 
respectively. By adding and substracting both equations 
(|l3|) , jl^) we obtain dynamical equations for two sets of 
dynamical hierarchies, each set characterized by a pop- 
ulation of oscillators with a given natural frequency (+1 
or -1). Note that the two sets of oscillators are also cou- 
pled one to each other through the terms ikg k and ikf k 
in (|l^),(|lj). These equations also can be used to ana- 
lytically compute stationary states and perform stability 
analysis as has been done in 113] . In this last case it 
can be shown that the fundamental model k = 1 decou- 
ples form the rest of the modes in a natural way. Here 
we follow a different strategy and use the method to nu- 
merically solve the set of equations for a given number 
2L + 1 of terms in the hierarchy {f k ,g k ; ~L < k < L). 
We stress that L is not the number of oscillators in the 
system which is already infinite from the beginning. To 
reduce possible dependences on the value of L we consider 
periodic boundary conditions f L+1 = fl L , f-L-i = ft 
and do the same for the g's. In our numerical solution 
we have taken L = 100 and we have checked that the 
results are the same by including more terms in the hier- 
archy. Equations have been solved using a second order 
Euler algorithm. We have also started from an initial 
condition of the type shown in eq. (|J) with 9 = in order 
to check they are not stationary solutions. Depending on 
the values of the parameters Kq and T there are different 
regimes | Q . 

In figures 2,3,4 we show the trajectory of the system 
in the plane (Re(fi), Im(g±)) for three different regimes 
(the incoherent, the critical and the coherent regimes). 
Note that according to eq.(fL2|) all the trajectories depart 
from the Im(gi) in the direction —i = exp(% c ). The 
first regime (fig-2) corresponds to the region where the 
incoherent solution is stable. In this case the order pa- 
rameter r (r — (/i/x)') oscillates with an amplitude 
which decays to zero exponentially in time. The second 
regime is shown in figure 3 and corresponds to the the 
critical boundary line T = ^ where the incoherent so- 
lution becomes unstable. In this case r oscillates and 
its amplitude decays to zero algebraically like t~ 2 as ex- 
pected for mean-field models at the critical point. In 
the region T < the incoherent solution is unstable 
and the system reaches a oscillating stationary solution 
(see figure 4) in agreement with the results analytically 
found by Bonilla et al fll3| , |l4|] . Note that this type of 
solutions cannot be computed from any theory based on 
an effective Hamiltonian (EH) like that given by eq.(|^). 
Therefore, the assumption of the existence of an effective 
Hamiltonian not only implies to consider solutions which 
are not stationary but also to miss another set of them 
that are explicitly time dependent. 
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FIG. 2. Trajectory of the system in the (Re(fi),Im(gi)) 
plane at T — 0.5, Kq = 1 in the incoherent regime. 
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FIG. 3. The same as in figure 2 at T = 0.25, K Q = 1 in the 
critical boundary. The area of the central hole decreases like 

We have also observed the existence of stationary fixed 
points for enough large values of Kq as expected when 
the ferromagnetic coupling is strong enough. In this case 
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the equilibrium solution within the EH approach albeit 
incorrect is closer to the true stationary one (in the limit 
K o — > oo the EH approach is recovered). 
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FIG. 4. The same as in figure 2 at T = 0.05, K = 1 in the 
synchronised regime. 

Finally, a comparison between the present theory 
(eqs. (|13yi4|) for the bimodal case) and the Brownian sim- 
ulations is shown in figure 5. Simulations have been per- 
formed by solving eq.(l) with an Euler method with a 
time step St = 0.005 and for a population of TV=50.000 
oscillators. Despite some small differences there is re- 
markable agreement between the simulation and analyt- 
ical results. 

In summary, we have presented an approach to the an- 
alytical solution of the Kuramoto model which is simple 
in its formulation and suitable for analytic and numerical 
computations. We have shown that the EH approach fails 
in predicting the stationary states of the system as well 
as the ground states of the energy function eq. ([?]). We 
can give a an explanation for this result. It has been sug- 
gested H that if a minimum of the energy function eq. (j?]) 
can be localized in the interior of the region [— tt, it) then 
this should be asymptotically stable. In the Kuramoto 
model it seems that this condition is indeed not satis- 
fied, at least for the ground states of eq.(Q). As shown 
in eq.(|l^) the ground states of eq.(0) are not stationary 
states of the dynamics. The quantitative violation of the 
stationarity property has been also analitically computed 
in eq.(|l2|). The reason for the discrepancy of our results 
with those reported in |(| relies on how the order of the 
limits N — > oo and t — > oo are taken. 
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FIG. 5. Analytical solution (lines) versus Brownian sim- 
ulations (points) in the oscillationg regime with parameters 
T = 2.5, Ko = 1/4. Brownian simulations were performed 
with 50000 oscillators and one realization of the noise. 

The limit t —> oo concerns particularly to the proper- 
ties of the stationary states. When the limit N — > oo 
is taken first then all ground states of Ti e f / become dy- 
namically unstable. The reason is that in this limit a 
finite density of oscillators in the ground state touch the 
boundaries where the effective Hamiltonian eq. (0) is dis- 
continuous. Obviously this is not true if the limit t — ► oo 
is taken first. Then, for finite N, the phases 0, of the 
oscillators do not touch the borders (because it has zero 
measure) and the ground states are stationary. As said 
before, the order of limits (first N — > oo and later t — > oo) 
is the one expected to describe the relevant asymptotic 
dynamics. We have used equations (0) to investigate the 
dynamical behavior of the Kuramoto model free of finite- 
size effects. Particular results have been obtained for the 
bimodal distribution but the method can be generally 
applied to any other distribution. It would be very in- 
teresting to use this approach to study the spectrum of 
correlation and response functions of the model as well as 
to investigate the presence of other dynamical regimes. 
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